# Install and load the necessary packages
if(!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse, rdrobust, readxl, viridis, grid, gridExtra, scales, cowplot, ggpubr, rddensity, rdd, sf, xtable, ggpattern, rmapshaper, kableExtra, ggh4x, gridGraphics)

# Load and prepare data ---------------------------------------------------

# Read the preprocessed data sets
districts_results <- read.csv("./data_coded/districts_results.csv")
survey_data <- read.csv("./data_coded/survey_data.csv")

# Select relevant respondents
representation_data <- survey_data %>% 
    # Select only respondents who voted with their Erststimme
  filter(erststimme > 0) %>% 
  
  # Remove respondents who voted for another party not belonging to the six major political parties in Germany 
  filter(erststimme != 801) %>% 
  
  # Create factor variables
  mutate(year = as.factor(year),
         party_name = as.factor(party_name)
         ) 

# Data set for main analysis: post-election respondents
post_data <- representation_data %>% filter(post == 1) %>% filter(complete.cases(rep_voters, rep_all))

# Data set for placebo analysis: pre-election respondents
pre_data <- representation_data %>% filter(post == 0) %>% filter(complete.cases(rep_voters, rep_all))

  


# Figure 1: Citizens' expectations to their local representative ----------

descriptive_own <- post_data %>% 
  # Make factor variable for better visualization
  mutate(rep_voters_fctr = factor(rep_voters,
                            levels = 1:5
  )) %>%
  
  # Total number of respondents
  mutate(n_total = n()) %>% 
  
  # Calculate proportions of each answer category
  group_by(rep_voters_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100) %>% 
  
  # Plot the results with a bar plot
  ggplot(aes(rep_voters_fctr, share)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", fill = viridis(3)[1]) +
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = 1:5,
    labels = c(
      "1 \nNot\nimportant\nat all",
      "2 \n\n",
      "3 \n\n",
      "4 \n\n",
      "5 \nVery\nimportant"
    )
  ) +
  scale_y_continuous(limits = c(0,56), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("Should represent\nown voters \nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )


# The code is the same as above; just for another outcome variable
descriptive_all <- post_data %>% 
  mutate(rep_all_fctr = factor(rep_all,
                            levels = 1:5
  )) %>%
  mutate(n_total = n()) %>% 
  group_by(rep_all_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100) %>% 
  ggplot(aes(rep_all_fctr, share)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", fill = viridis(3)[2]) +
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = 1:5,
    labels = c(
      "1 \nNot\nimportant\nat all",
      "2 \n\n",
      "3 \n\n",
      "4 \n\n",
      "5 \nVery\nimportant"
    )
  ) +
  scale_y_continuous(limits = c(0,56), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("Should represent\nall citizens \nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )



# The code is the same as above; just for another outcome variable
descriptive_relative <- post_data %>% 
  mutate(relative_importance_fctr = factor(relative_importance,
                                           levels = -4:4
  )) %>%
  mutate(n_total = n()) %>% 
  group_by(relative_importance_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100) %>% 
  
  mutate(above0 = case_when(
    relative_importance_fctr %in% c(-4,-3,-2,-1,0) ~ "no",
    relative_importance_fctr %in% c(1,2,3,4) ~ "yes"
  )) %>% 
  
  ggplot(aes(relative_importance_fctr, share)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", aes(fill = above0)) +
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = -4:4
  ) +
  scale_y_continuous(limits = c(0,56), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("\nRelative\nimportance") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  scale_fill_manual("", values = c("no" = viridis(3)[2], "yes" = viridis(3)[1])) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )


ggarrange(descriptive_own, NULL, descriptive_all, NULL, descriptive_relative,
          ncol = 5, nrow = 1, align = "hv",
          widths = c(2,0.05,2,0.05,2)
)

ggsave("plots/figure_1.pdf", width = 2100, height = 900, units = "px")







# Figure 2: Effect of winning/losing on representative conceptions (regression discontinuity design) --------------------------------------------------------------

# Set names of different model specifications 
model_specification <- c("Conventional", "Bias corrected", "Robust")

# RDD estimation: Represent own voters 
rdd_own <- rdrobust::rdrobust(y = post_data$rep_voters, 
                              x = post_data$margin, 
                              c = 0, 
                              p = 1,
                              q = 2,
                              covs = cbind(post_data$party_name, post_data$year), 
                              vce = "nn", 
                              cluster = post_data$cluster, 
                              all = TRUE
)

# Create empty data frame in which the results will be piped in
results_own <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)

# Pipe in the results from the RDD models
for(i in 1: length(model_specification)){
  
  results_own[i, "specification"] <- model_specification[i]
  results_own[i, "LATE"] <- rdd_own$coef[i] 
  results_own[i, "lower"] <- rdd_own$ci[i,1]
  results_own[i, "upper"] <- rdd_own$ci[i,2]
  results_own[i, "SE"] <- rdd_own$se[i]
  results_own[i, "p"] <- rdd_own$pv[i]
  results_own[i, "bw_est"] <- rdd_own$bws[1,1]
  results_own[i, "bw_bias"] <- rdd_own$bws[2,1]
  results_own[i, "total_obs"] <- rdd_own$N %>% sum()
  results_own[i, "eff_obs"] <- if_else(i == 1, rdd_own$N_h %>% sum(), rdd_own$N_b %>% sum())
  results_own[i, "left"] <- if_else(i == 1, rdd_own$N_h[1], rdd_own$N_b[1])
  results_own[i, "right"] <- if_else(i == 1, rdd_own$N_h[2], rdd_own$N_b[2])
}


# Change labels of model specifciation for plotting
results_own <- results_own %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))



# Note: the code below here repeats the same steps, just for different outcome variables

# RDD estimation: Represent all citizens 
rdd_all <- rdrobust::rdrobust(y = post_data$rep_all, 
                              x = post_data$margin, 
                              c = 0, 
                              p = 1,
                              q = 2,
                              covs = cbind(post_data$party_name, post_data$year), 
                              vce = "nn", 
                              cluster = post_data$cluster, 
                              all = TRUE
)


model_specification <- c("Conventional", "Bias corrected", "Robust")

results_all <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_all[i, "specification"] <- model_specification[i]
  results_all[i, "LATE"] <- rdd_all$coef[i] 
  results_all[i, "lower"] <- rdd_all$ci[i,1]
  results_all[i, "upper"] <- rdd_all$ci[i,2]
  results_all[i, "SE"] <- rdd_all$se[i]
  results_all[i, "p"] <- rdd_all$pv[i]
  results_all[i, "bw_est"] <- rdd_all$bws[1,1]
  results_all[i, "bw_bias"] <- rdd_all$bws[2,1]
  results_all[i, "total_obs"] <- rdd_all$N %>% sum()
  results_all[i, "eff_obs"] <- if_else(i == 1, rdd_all$N_h %>% sum(), rdd_all$N_b %>% sum())
  results_all[i, "left"] <- if_else(i == 1, rdd_all$N_h[1], rdd_all$N_b[1])
  results_all[i, "right"] <- if_else(i == 1, rdd_all$N_h[2], rdd_all$N_b[2])
}


results_all <- results_all %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))



# RDD estimation: Relative importance
rdd_relative <- rdrobust::rdrobust(y = post_data$relative_importance, 
                                   x = post_data$margin, 
                                   c = 0, 
                                   p = 1,
                                   q = 2,
                                   covs = cbind(post_data$party_name, post_data$year), 
                                   vce = "nn", 
                                   cluster = post_data$cluster, 
                                   all = TRUE
)

model_specification <- c("Conventional", "Bias corrected", "Robust")

results_relative <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_relative[i, "specification"] <- model_specification[i]
  results_relative[i, "LATE"] <- rdd_relative$coef[i] 
  results_relative[i, "lower"] <- rdd_relative$ci[i,1]
  results_relative[i, "upper"] <- rdd_relative$ci[i,2]
  results_relative[i, "SE"] <- rdd_relative$se[i]
  results_relative[i, "p"] <- rdd_relative$pv[i]
  results_relative[i, "bw_est"] <- rdd_relative$bws[1,1]
  results_relative[i, "bw_bias"] <- rdd_relative$bws[2,1]
  results_relative[i, "total_obs"] <- rdd_relative$N %>% sum()
  results_relative[i, "eff_obs"] <- if_else(i == 1, rdd_relative$N_h %>% sum(), rdd_relative$N_b %>% sum())
  results_relative[i, "left"] <- if_else(i == 1, rdd_relative$N_h[1], rdd_relative$N_b[1])
  results_relative[i, "right"] <- if_else(i == 1, rdd_relative$N_h[2], rdd_relative$N_b[2])
}



results_relative <- results_relative %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))






### Own voters
# Plot the coefficients
rdd_own_coefs <- results_own %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(-0.3,0.5) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11)) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )


# RDD plot:

# Set minimum and maximum values for plotting on the y-axis scale
ymin <- 3
ymax <- 5

# Create an RDD plot
rdd_own_rdplot <- rdplot(
  y = post_data$rep_voters,
  x = post_data$margin,
  c = 0,
  covs = cbind(post_data$party_name, post_data$year),
  #weights = post_data$weight,
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  p = 2,
  title = "",
  x.label = "Margin",
  y.label = "Importance"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin, ymax), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("Should represent\nown voters\nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_own$bw_est[1],
    xmax = results_own$bw_est[1],
    ymin = ymin,
    ymax = ymax
  ),
  fill = "grey",
  alpha = 0.5)

### All voters
# Plot the coefficients
rdd_all_coefs <- results_all %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(-0.3,0.5) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11))+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )
# RDD plot:

# Set minimum and maximum values for plotting on the y-axis scale
ymin <- 3
ymax <- 5

# Create an RDD plot
rdd_all_rdplot <- rdplot(
  y = post_data$rep_all,
  x = post_data$margin,
  c = 0,
  covs = cbind(post_data$party_name, post_data$year),
  #weights = post_data$weight,
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  p = 2,
  title = "",
  x.label = "Margin",
  y.label = "Importance"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin, ymax), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("Should represent\nall citizens \nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_all$bw_est[1],
    xmax = results_all$bw_est[1],
    ymin = ymin,
    ymax = ymax
  ),
  fill = "grey",
  alpha = 0.5)

### Relative importance
# Plot the coefficients
rdd_relative_coefs <- results_relative %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(-0.3,0.5) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11))+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )
# RDD plot:

# Set minimum and maximum values for plotting on the y-axis scale
ymin_relative <- -0.5
ymax_relative <- 1

# Create an RDD plot
rdd_relative_rdplot <- rdplot(
  y = post_data$relative_importance,
  x = post_data$margin,
  c = 0,
  covs = cbind(post_data$party_name, post_data$year),
  #weights = post_data$weight,
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  p = 2,
  title = "",
  x.label = "Margin",
  y.label = "Relative importance"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin_relative, ymax_relative), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("\nRelative \nimportance") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_relative$bw_est[1],
    xmax = results_relative$bw_est[1],
    ymin = ymin_relative,
    ymax = ymax_relative
  ),
  fill = "grey",
  alpha = 0.5)


# Combine all plots together
plot_grid(
  rdd_own_rdplot,
  rdd_all_rdplot,
  rdd_relative_rdplot,
  rdd_own_coefs,
  rdd_all_coefs,
  rdd_relative_coefs,
  
  align = "v",
  nrow = 2 ,   
  rel_heights = c(2.3/3, 1/3)
)

ggsave("plots/figure_2.pdf", width = 2100, height = 1250, units = "px")







# Online Appendix ---------------------------------------------------------

## Table C1: RDD full results of Figure 2  ---------------------------------

# Create a table displying the results from the three main RDD models
main <- results_own %>% 
  mutate(dv = "Own voters") %>% 
  mutate(specification = factor(specification,
                                c("Conven-\ntional", "Bias \ncorrec.", "Robust"),
                                labels = c("Conventional", "Bias corrected", "Robust")
  )) %>% 
  select(dv, specification, LATE, SE, p, bw_est, bw_bias, total_obs, eff_obs, left, right) %>% 
  rbind(
    results_all %>% 
      mutate(dv = "All constituents") %>% 
      mutate(specification = factor(specification,
                                    c("Conven-\ntional", "Bias \ncorrec.", "Robust"),
                                    labels = c("Conventional", "Bias corrected", "Robust")
      )) %>% 
      select(dv, specification, LATE, SE, p, bw_est, bw_bias, total_obs, eff_obs, left, right) %>% 
      rbind(
        results_relative %>% 
          mutate(dv = "Relative importance") %>% 
          mutate(specification = factor(specification,
                                        c("Conven-\ntional", "Bias \ncorrec.", "Robust"),
                                        labels = c("Conventional", "Bias corrected", "Robust")
          )) %>% 
          select(dv, specification, LATE, SE, p, bw_est, bw_bias, total_obs, eff_obs, left, right)
      )
  ) %>% 
  rename(
    bw.est = bw_est,
    bw.bias = bw_bias,
    total.obs = total_obs,
    eff.obs = eff_obs
  ) %>% 
  dplyr::mutate_at(
    .vars = vars(c("LATE", "SE", "p")),
    .funs = ~ format(round(., digits = 3), nsmall = 3)
  ) %>% 
  mutate(bw = dplyr::if_else(specification == "Conventional", bw.est, bw.bias)) %>% 
  mutate(bw = format(round(bw, digits = 2), nsmall = 2)) %>% 
  select(specification, LATE, SE, p, bw, total.obs, eff.obs, left, right)

# Transpose matrix
main_table <- t(main)

# Change column names 
colnames(main_table) <- main_table[1,]

main_table <- main_table[-1,]

rownames(main_table) <- c("LATE", "SE", "P>|z|", "Bandwidth", "Total no. obs.", "Effect. no. obs.", "Left of cutoff", "Right of cutoff")

# Print out the table so it can be put into latex
kableExtra::kable(
  main_table,
  format = "latex",
  booktabs = TRUE,
  escape = FALSE,
  align = "c",
  caption = "RDD full results of Figure \\ref{fig:main-plot} in the main part.",
  linesep = ""
)  %>% 
  kableExtra::kable_styling(font_size = 10, full_width = TRUE) %>% 
  kableExtra::row_spec(3, hline_after = TRUE) %>% 
  kableExtra::add_header_above(header = c(" " = 1, "Own voters in constituency" = 3, "All citizens in constituency" = 3, "Relative importance" = 3), escape = FALSE) %>% 
  kableExtra::landscape()





## Figure A1.1: Winning candidates in electoral districts (Erststimme) --------

# Read the shape files of the electoral districts in each year
wahlkreise_shp_13 <- read_sf("raw_data/btw13_geometrie_wahlkreise_shp/Geometrie_Wahlkreise_18DBT.shp", quiet = TRUE) 
wahlkreise_shp_17 <- read_sf("raw_data/btw17_geometrie_wahlkreise_shp/Geometrie_Wahlkreise_19DBT.shp", quiet = TRUE) 
wahlkreise_shp_21 <- read_sf("raw_data/btw21_geometrie_wahlkreise_shp/Geometrie_Wahlkreise_20DBT.shp", quiet = TRUE) 


district_winners <- districts_results %>% 
  
  # Only get rows of winning party in each district and election 
  filter(party_name == winner_name) %>% 
  
  # Create variable indicating whether district race was close (margin of victory smaller than 5 percentage points) for plotting
  mutate(close = dplyr::if_else(margin <= 5, "close", "clear")) %>%
  
  # Select relevant variables
  select(year, wahlkreis, winner_name, margin, close) %>% 
  
  # Identify how many district races each party won in each year
  group_by(winner_name, year) %>% 
  mutate(n_won = n()) %>% 
  ungroup() %>% 
  
  # Create new party labels for plotting (with info how many district races they won)
  mutate(winner_name_full = factor(winner_name,
                                   levels = c("cdu", "spd","fdp", "green", "left", "afd", "others"),
                                   labels = c("CDU/CSU", "SPD", "FDP", "Greens", "The Left", "AfD", "Others") 
  ))  %>% 
  
  mutate(winner_label = paste(winner_name_full, " (N = ", n_won, ")", sep = "")) 


# Map of district winners: 2013
map_2013 <- wahlkreise_shp_13 %>% 
  
  # Filter 2013 years and set factor levels of the labels
  left_join(., district_winners %>% filter(year == 2013) %>% 
              
              mutate(winner_label = factor(winner_label, levels = c("CDU/CSU (N = 236)",
                                                                    "SPD (N = 58)",
                                                                    "Greens (N = 1)",
                                                                    "The Left (N = 4)"
              )) )
            
            , by = c("WKR_NR" = "wahlkreis")) %>% 
  
  # Simplify polygons slightly
  rmapshaper::ms_simplify(., keep = 0.1,
                          keep_shapes = FALSE) %>%
  
  # PLot the map
  ggplot() +
  
  geom_sf_pattern(
    aes(fill = winner_label, pattern = close),
    pattern_colour = "white",   # Color of the pattern
    pattern_fill="white",
    pattern_density = 0.01,    # Adjust pattern density
    pattern_spacing = 0.01,   # Spacing of stripes
    lwd = 0.03,
    show.legend = c(pattern = FALSE)
  ) +
  
  # Pattern indicates whether district race was close or not
  scale_pattern_manual(values = c("none", "pch"), guide = "none") +
  
  # Colors for the parties
  scale_fill_manual(values = c("black", "red", "green", "purple"), name = NULL,
                    
                    guide_legend(override.aes = list(pattern = 'none'))
  ) +
  
  theme_bw() +
  coord_sf(datum = NA) +
  theme(panel.border=element_blank(),
        panel.grid = element_blank(),
        axis.text.x= element_blank(),
        axis.text.y = element_blank()) +
  
  ggtitle("2013") +
  theme(plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
        legend.text = element_text(size = 15),
        legend.position = "bottom"
  ) +
  guides(fill=guide_legend(nrow=5,byrow=TRUE, override.aes = list(pattern = "none"))
  ) 


# NOTE: The code below follows the same structure as above; just for different years

# Map of district winners: 2017
map_2017 <- wahlkreise_shp_17 %>% 
  left_join(., district_winners %>% filter(year == 2017) %>% 
              
              mutate(winner_label = factor(winner_label, levels = c("CDU/CSU (N = 231)",
                                                                    "SPD (N = 59)",
                                                                    "Greens (N = 1)",
                                                                    "The Left (N = 5)",
                                                                    "AfD (N = 3)"
              )) )
            
            , by = c("WKR_NR" = "wahlkreis")) %>% 
  
  # Simplify polygons slightly
  rmapshaper::ms_simplify(., keep = 0.1,
                          keep_shapes = FALSE) %>%
  
  ggplot() +
  
  geom_sf_pattern(
    aes(fill = winner_label, pattern = close),
    pattern_colour = "white",   # Color of the pattern
    pattern_fill="white",
    pattern_density = 0.01,    # Adjust pattern density
    pattern_spacing = 0.01,   # Spacing of stripes
    lwd = 0.03,
    show.legend = c(pattern = FALSE)
  ) +
  
  scale_pattern_manual(values = c("none", "pch"), guide = "none") +
  
  
  scale_fill_manual(values = c("black", "red", "green", "purple", "deepskyblue"), name = NULL) +
  theme_bw() +
  coord_sf(datum = NA) +
  theme(panel.border=element_blank(),
        panel.grid = element_blank(),
        axis.text.x= element_blank(),
        axis.text.y = element_blank()) +
  
  ggtitle("2017") +
  theme(plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
        legend.text = element_text(size = 15),
        legend.position = "bottom"
  )+
  guides(fill=guide_legend(nrow=5,byrow=TRUE, override.aes = list(pattern = "none"))
  ) 






# Map of district winners: 2021
map_2021 <- wahlkreise_shp_21 %>% 
  left_join(., district_winners %>% filter(year == 2021) %>% 
              
              mutate(winner_label = factor(winner_label, levels = c("CDU/CSU (N = 143)",
                                                                    "SPD (N = 121)",
                                                                    "Greens (N = 16)",
                                                                    "The Left (N = 3)",
                                                                    "AfD (N = 16)"
              )) )
            
            , by = c("WKR_NR" = "wahlkreis")) %>% 
  
  # Simplify polygons slightly
  rmapshaper::ms_simplify(., keep = 0.1,
                          keep_shapes = FALSE) %>%
  
  ggplot() +
  
  
  geom_sf_pattern(
    aes(fill = winner_label, pattern = close),
    pattern_colour = "white",   # Color of the pattern
    pattern_fill="white",
    pattern_density = 0.01,    # Adjust pattern density
    pattern_spacing = 0.01,   # Spacing of stripes
    lwd = 0.03,
    show.legend = c(pattern = FALSE)
  ) +
  
  scale_pattern_manual(values = c("none", "pch"), guide = "none") +
  
  scale_fill_manual(values = c("black", "red", "green", "purple", "deepskyblue"), name = NULL) +
  theme_bw() +
  coord_sf(datum = NA) +
  theme(panel.border=element_blank(),
        panel.grid = element_blank(),
        axis.text.x= element_blank(),
        axis.text.y = element_blank()) +
  
  ggtitle("2021") +
  theme(plot.title = element_text(size = 20, face = "bold", hjust = 0.5),
        legend.text = element_text(size = 15),
        legend.position = "bottom"
  ) +
  guides(fill=guide_legend(nrow=5,byrow=TRUE, override.aes = list(pattern = "none"))
  ) 



# Combine all the maps together
map_all <- cowplot::plot_grid(map_2013, NULL, map_2017, NULL, map_2021, 
                              nrow = 1,
                              rel_widths = c(1,0.1,1,0.1,1)
)


ggsave(map_all, filename = "plots/figure_a1.1.pdf", width = 12, height = 8)





## Figure A2.1: Election results in districts ------------------------------


margin_plot <- districts_results %>% 
  filter(party == 1 & year != 2009) %>% 
  mutate(margin = winner_vote_share - second_vote_share) %>% 
  mutate(close = if_else(margin < 5, 1, 0)) %>% 
  ggplot(aes(x = margin)) +
  theme_bw() +
  geom_histogram(binwidth = 0.5, fill = viridis(1)[1], color = viridis(1)[1], size = 0.2) +
  scale_x_continuous(limits = c(-1,52), expand = c(0,0)) +  
  scale_y_continuous(limits = c(0,28), expand = c(0,0), breaks = seq(5,25,5)) +
  # ggtitle("Overview of the closeness of the district elections") +
  # theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Closeness \n(in percentage points)") +
  ylab("Number of elections") +
  ggtitle("Closeness of the \ndistrict elections") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  )



winner_share_plot <- districts_results %>% 
  filter(party == 1 & year != 2009) %>% 
  mutate(margin_winner = winner_vote_share - second_vote_share) %>% 
  filter(margin_winner <= 2.5) %>% 
  ggplot(aes(x = winner_vote_share)) +
  theme_bw() +
  geom_histogram(binwidth = 1, fill = viridis(1)[1], color = viridis(1)[1], size = 0.2) +
  scale_x_continuous(limits = c(-1,52), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,13), expand = c(0,0), breaks = seq(2,12,2)) +
  # ggtitle("Overview of the closeness of the district elections") +
  # theme(plot.title = element_text(hjust = 0.5)) +
  xlab("Vote share of winning candidate \n(in percent)") +
  ylab("Number of elections") +
  ggtitle("Vote share of winning \ncandidates in close elections") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  )

cowplot::plot_grid(margin_plot, winner_share_plot)

ggsave("plots/figure_a2.1.pdf", width = 1800, height = 900, units = "px")





## Figure A2.2: Margin of victory/defeat in district races by political party affiliaition of candidates --------


districts_results %>% 
  
  # Set factor levels of parties for plotting
  mutate(party_name = factor(party_name,
                             levels = c("cdu", "spd", "fdp", "green", "left", "afd"),
                             labels = c("CDU/CSU", "SPD", "FDP", "Greens", "The Left", "AfD")
  )) %>% 
  
  # Plot histograms for each party
  ggplot() +
  
  geom_histogram(aes(x = margin), , fill = viridis(1)[1], color = viridis(1)[1], size = 0.2, center = 0.5, binwidth = 1) +
  
  facet_wrap(~ party_name, nrow = 2) +
  
  theme_bw() +
  
  scale_y_continuous(limits = c(0, 43), expand = c(0,0)) +
  
  geom_vline(xintercept = 0,
             lty = 2,
             linewidth = 0.25,
             colour = gray(1 / 2)) +
  
  xlab("Margin \n(in percentage points)") +
  ylab("Number of elections") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 14),
        strip.text = element_text(face = "bold", size = 12),
        axis.text = element_text(size = 11),
        axis.title = element_text(size = 12)
  ) +
  theme(panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  )

ggsave("plots/figure_a2.2.pdf", width = 2650, height = 1700, units = "px")





## Table B1: Distribution of perceived importances -------------------------

missing_data <- representation_data %>%
  
  # Filter post election responses (need to build up on other data.frame, as missing responses are filtered out in "post_data" object)
  filter(post == 1) %>% 
  
  # Filter out those who stopped the interview early (didn't reach these question items)
  filter(rep_voters_raw!= -93) %>% 
  
  # New factor variable: now also includes values that are coded as missing in the main analysis
  mutate(rep_voters_original = case_when(
    rep_voters_raw == 1 ~ "Very important",
    rep_voters_raw == 2 ~ "Somewhat important",
    rep_voters_raw == 3 ~ "In between",
    rep_voters_raw == 4 ~ "Not very important",
    rep_voters_raw == 5 ~ "Not important at all",
    rep_voters_raw == -98 ~ "Don't know",
    rep_voters_raw == -99 ~ "No answer"
  )) %>% 
  
  mutate(rep_voters_original = factor(rep_voters_original,
                                levels = c("Very important",
                                           "Somewhat important",
                                           "In between",
                                           "Not very important",
                                           "Not important at all",
                                           "Don't know",
                                           "No answer")
  )) %>% 
  
  # Here: same as above for other variable
  mutate(rep_all_original = case_when(
    rep_all_raw == 1 ~ "Very important",
    rep_all_raw == 2 ~ "Somewhat important",
    rep_all_raw == 3 ~ "In between",
    rep_all_raw == 4 ~ "Not very important",
    rep_all_raw == 5 ~ "Not important at all",
    rep_all_raw == -98 ~ "Don't know",
    rep_all_raw == -99 ~ "No answer"
  )) %>% 
  
  mutate(rep_all_original = factor(rep_all_original,
                                levels = c("Very important",
                                           "Somewhat important",
                                           "In between",
                                           "Not very important",
                                           "Not important at all",
                                           "Don't know",
                                           "No answer")
  )) 


# Table of responses for outcome: Should represent own voters in constituency
rep_voters_counts <- missing_data$rep_voters_original %>% 
  table() %>% 
  as.data.frame() %>% 
  rename(category = ".",
         n = Freq
  ) %>% 
  mutate(n_total = sum(n),
         percent = n/n_total * 100
  ) %>% 
  
  mutate(out = paste(n, " (", round(percent, digits = 2), "%)", sep = "")) %>% 
  select(category, out) %>% 
  
  pivot_wider(
    names_from = category,
    values_from = out
  ) %>% 
  
  mutate(Item = "Should represent own voters in constituency") %>% 
  select(Item, `Very important`, `Somewhat important`, `In between`, `Not very important`, `Not important at all`, `Don't know`, `No answer`)


# Table of responses for outcome: Should represent all citizens in constituency
rep_all_counts <- missing_data$rep_all_original %>% 
  table() %>% 
  as.data.frame() %>% 
  rename(category = ".",
         n = Freq
  ) %>% 
  mutate(n_total = sum(n),
         percent = n/n_total * 100
  ) %>% 
  
  mutate(out = paste(n, " (", round(percent, digits = 2), "%)", sep = "")) %>% 
  select(category, out) %>% 
  
  pivot_wider(
    names_from = category,
    values_from = out
  ) %>% 
  
  mutate(Item = "Should represent all citizens in constituency") %>% 
  select(Item, `Very important`, `Somewhat important`, `In between`, `Not very important`, `Not important at all`, `Don't know`, `No answer`)


print(xtable::xtable(rbind(rep_voters_counts, rep_all_counts)), include.rownames=FALSE)







## Figure B2.1: Perceived importance of different representational foci --------

# Bar plot: Should represent own voters in constituency
descriptive_own <- post_data %>% 
  filter(!is.na(rep_voters)) %>% 
  mutate(rep_voters_fctr = factor(rep_voters,
                            levels = 1:5
  )) %>%
  mutate(n_total = n()) %>% 
  group_by(rep_voters_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100,
         share_round = paste(formatC(round(share, digits = 2), format = "f", digits = 2), "%", sep = "")
  ) %>% 
  ggplot(aes(rep_voters_fctr, share, label = share_round)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", fill = viridis(3)[1]) +
  
  geom_text(position = position_dodge(width = .9),    # move to center of bars
            vjust = -0.5,    # nudge above top of bar
            size = 3) +
  
  
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = 1:5,
    labels = c(
      "1 \nNot\nimportant\nat all",
      "2 \n\n",
      "3 \n\n",
      "4 \n\n",
      "5 \nVery\nimportant"
    )
  ) +
  scale_y_continuous(limits = c(0,61), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("Should represent own\nvoters in constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )


# Bar plot: Should represent all citizens in constituency
descriptive_all <- post_data %>% 
  filter(!is.na(rep_all)) %>% 
  mutate(rep_all_fctr = factor(rep_all,
                            levels = 1:5
  )) %>%
  mutate(n_total = n()) %>% 
  group_by(rep_all_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100,
         share_round = paste(formatC(round(share, digits = 2), format = "f", digits = 2), "%", sep = "")
  ) %>% 
  ggplot(aes(rep_all_fctr, share, label = share_round)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", fill = viridis(3)[1]) +
  
  geom_text(position = position_dodge(width = .9),    # move to center of bars
            vjust = -0.5,    # nudge above top of bar
            size = 3) +
  
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = 1:5,
    labels = c(
      "1 \nNot\nimportant\nat all",
      "2 \n\n",
      "3 \n\n",
      "4 \n\n",
      "5 \nVery\nimportant"
    )
  ) +
  scale_y_continuous(limits = c(0,61), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("Should represent all\ncitizens in constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )



# Bar plot: Should represent all voters of own party
descriptive_party <- 
  post_data %>% as.data.frame() %>% 
  filter(!is.na(rep_party)) %>% 
  mutate(rep_party_fctr = factor(rep_party,
                            levels = 1:5
  )) %>%
  mutate(n_total = n()) %>% 
  group_by(rep_party_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100,
         share_round = paste(formatC(round(share, digits = 2), format = "f", digits = 2), "%", sep = "")
  ) %>% 
  ggplot(aes(rep_party_fctr, share, label = share_round)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", fill = viridis(3)[1]) +
  
  geom_text(position = position_dodge(width = .9),    # move to center of bars
            vjust = -0.5,    # nudge above top of bar
            size = 3) +
  
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = 1:5,
    labels = c(
      "1 \nNot\nimportant\nat all",
      "2 \n\n",
      "3 \n\n",
      "4 \n\n",
      "5 \nVery\nimportant"
    )
  ) +
  scale_y_continuous(limits = c(0,61), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("Should represent all\nvoters of own party") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  ) +
  theme(plot.margin = unit(c(0.1,0.3,-0.1,0.3), "cm"))




# Bar plot: Should represent all citizens in Germany
descriptive_nation <- 
  post_data %>% as.data.frame() %>% 
  filter(!is.na(rep_nation)) %>% 
  mutate(rep_nation_fctr = factor(rep_nation,
                            levels = 1:5
  )) %>%
  mutate(n_total = n()) %>% 
  group_by(rep_nation_fctr, n_total) %>% 
  summarise(n = n()) %>% 
  rowwise() %>% 
  mutate(share = n / n_total * 100,
         share_round = paste(formatC(round(share, digits = 2), format = "f", digits = 2), "%", sep = "")
  ) %>% 
  ggplot(aes(rep_nation_fctr, share, label = share_round)) + 
  theme_bw() +
  geom_bar(position = "dodge", stat="identity", fill = viridis(3)[1]) +
  
  geom_text(position = position_dodge(width = .9),    # move to center of bars
            vjust = -0.5,    # nudge above top of bar
            size = 3) +
  
  ylab("% of respondents") +
  xlab(element_blank()) +
  scale_x_discrete(
    breaks = 1:5,
    labels = c(
      "1 \nNot\nimportant\nat all",
      "2 \n\n",
      "3 \n\n",
      "4 \n\n",
      "5 \nVery\nimportant"
    )
  ) +
  scale_y_continuous(limits = c(0,61), expand = c(0,0)) +
  theme(legend.position="none") +
  theme(legend.title=element_blank()) +
  ggtitle("Should represent all\ncitizens in Germany") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10),
        axis.title = element_text(size = 11)
  ) +
  
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  ) +
  theme(plot.margin = unit(c(0.1,0.3,-0.1,0.3), "cm"))


# Combine plots together
cowplot::plot_grid(descriptive_own,
                   descriptive_all,
                   descriptive_party,
                   descriptive_nation,
                   nrow = 2
)



ggsave("plots/figure_b2.1.pdf", width = 2300, height = 1600, units = "px")






## Table B3: Vote choices of respondents -----------------------------------


# Use raw data set here: because don'T knows etc. got filtered out for the main analysis 
votechoice_data <- survey_data %>%
  
  # Select post-election respondents
  filter(post == 1) %>% 
  # Filter out respondents for which this question did not apply (e.g. due to no turnout) or who did not reach this question
  filter(!(erststimme_raw %in% c(-97, -73))) %>% 
  filter(!(zweitstimme_raw %in% c(-97, -73))) %>% 
  
  # Labels for votes
  mutate(erststimme_full = dplyr::recode(erststimme_raw,
                                         "1" = "CDU/CSU",
                                         "4" = "SPD",
                                         "5" = "FDP",
                                         "6" = "Greens",
                                         "7" = "The Left",
                                         "322" = "AfD",
                                         "801" = "Others",
                                         "206" = "Others",
                                         "215" = "Others",
                                         "-99" = "No answer",
                                         "-98" = "Don't know",
                                         "-83" = "Invalid vote"
  )) %>% 
  mutate(zweitstimme_full = dplyr::recode(zweitstimme_raw,
                                          "1" = "CDU/CSU",
                                          "4" = "SPD",
                                          "5" = "FDP",
                                          "6" = "Greens",
                                          "7" = "The Left",
                                          "322" = "AfD",
                                          "801" = "Others",
                                          "206" = "Others",
                                          "215" = "Others",
                                          "-99" = "No answer",
                                          "-98" = "Don't know",
                                          "-83" = "Invalid vote"
  )) %>% 
  
  select(year, erststimme_full, zweitstimme_full)




# Create table displaying reported vote choices of respondents
vote_choices <- votechoice_data %>% 
  
  pivot_longer(cols = c("erststimme_full", "zweitstimme_full"),
               names_to = "vote_type",
               values_to = "party") %>% 
  
  group_by(year, vote_type) %>% 
  
  mutate(n_total = n()) %>% 
  ungroup() %>% 
  
  
  mutate(party = factor(party,
                        levels = c(
                          "CDU/CSU",
                          "SPD",
                          "FDP",
                          "Greens",
                          "The Left",
                          "AfD",
                          "Others",
                          "Invalid vote",
                          "Don't know",
                          "No answer"
                        )
  )) %>% 
  
  group_by(year, vote_type, party, n_total) %>% 
  count() %>% 
  ungroup() %>% 
  
  pivot_wider(names_from = vote_type, values_from = n) %>% 
  
  mutate(erststimme_percent = erststimme_full / n_total * 100,
         zweitstimme_percent = zweitstimme_full / n_total * 100
  ) %>% 
  
  mutate(out_erststimme = paste(erststimme_full, " (", formatC(round(erststimme_percent, digits = 2), format = "f", digits = 2), "%)", sep = "")) %>% 
  
  mutate(out_zweitstimme = paste(zweitstimme_full, " (", formatC(round(zweitstimme_percent, digits = 2), format = "f", digits = 2), "%)", sep = "")) %>% 
  
  select(year, party, out_erststimme, out_zweitstimme) %>% 
  
  pivot_wider(names_from = year, values_from = c(out_erststimme, out_zweitstimme)) %>% 
  select(party, 
         out_erststimme_2013, 
         out_zweitstimme_2013, 
         out_erststimme_2017, 
         out_zweitstimme_2017,
         out_erststimme_2021,
         out_zweitstimme_2021
  ) %>% 
  
  rename(Party = party,
         `First 2013` = out_erststimme_2013,
         `Second 2013` = out_zweitstimme_2013, 
         `First 2017` = out_erststimme_2017, 
         `Second 2017` = out_zweitstimme_2017,
         `First 2021` = out_erststimme_2021,
         `Second 2021` = out_zweitstimme_2021
  )



print(xtable::xtable(vote_choices), include.rownames=FALSE)







## Figure D1.1: Effect of winning/losing on alternative operationalizations of the relative importance measure -------------------------------------------------

# NOTE: The code follows the exact same structure as the code for Figure 2; just with different outcome variables

# Set names of different model specification
model_specification <- c("Conventional", "Bias corrected", "Robust")

# RDD: Binary outcome measure
rdd_binary <- rdrobust::rdrobust(y = post_data$relative_binary, 
                                 x = post_data$margin, 
                                 c = 0, 
                                 p = 1,
                                 q = 2,
                                 #h = 1,
                                 covs = cbind(post_data$party_name, post_data$year), 
                                 vce = "nn", 
                                 cluster = post_data$cluster, 
                                 all = TRUE
)

model_specification <- c("Conventional", "Bias corrected", "Robust")

results_binary <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_binary[i, "specification"] <- model_specification[i]
  results_binary[i, "LATE"] <- rdd_binary$coef[i] 
  results_binary[i, "lower"] <- rdd_binary$ci[i,1]
  results_binary[i, "upper"] <- rdd_binary$ci[i,2]
  results_binary[i, "SE"] <- rdd_binary$se[i]
  results_binary[i, "p"] <- rdd_binary$pv[i]
  results_binary[i, "bw_est"] <- rdd_binary$bws[1,1]
  results_binary[i, "bw_bias"] <- rdd_binary$bws[2,1]
  results_binary[i, "total_obs"] <- rdd_binary$N %>% sum()
  results_binary[i, "eff_obs"] <- if_else(i == 1, rdd_binary$N_h %>% sum(), rdd_binary$N_b %>% sum())
  results_binary[i, "left"] <- if_else(i == 1, rdd_binary$N_h[1], rdd_binary$N_b[1])
  results_binary[i, "right"] <- if_else(i == 1, rdd_binary$N_h[2], rdd_binary$N_b[2])
}


results_binary <- results_binary %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))



# RDD: Ratio outcome measure
rdd_ratio <- rdrobust::rdrobust(y = post_data$relative_ratio, 
                                x = post_data$margin, 
                                c = 0, 
                                p = 1,
                                q = 2,
                                #h = 1,
                                covs = cbind(post_data$party_name, post_data$year), 
                                vce = "nn", 
                                cluster = post_data$cluster, 
                                all = TRUE
)

results_ratio <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_ratio[i, "specification"] <- model_specification[i]
  results_ratio[i, "LATE"] <- rdd_ratio$coef[i] 
  results_ratio[i, "lower"] <- rdd_ratio$ci[i,1]
  results_ratio[i, "upper"] <- rdd_ratio$ci[i,2]
  results_ratio[i, "SE"] <- rdd_ratio$se[i]
  results_ratio[i, "p"] <- rdd_ratio$pv[i]
  results_ratio[i, "bw_est"] <- rdd_ratio$bws[1,1]
  results_ratio[i, "bw_bias"] <- rdd_ratio$bws[2,1]
  results_ratio[i, "total_obs"] <- rdd_ratio$N %>% sum()
  results_ratio[i, "eff_obs"] <- if_else(i == 1, rdd_ratio$N_h %>% sum(), rdd_ratio$N_b %>% sum())
  results_ratio[i, "left"] <- if_else(i == 1, rdd_ratio$N_h[1], rdd_ratio$N_b[1])
  results_ratio[i, "right"] <- if_else(i == 1, rdd_ratio$N_h[2], rdd_ratio$N_b[2])
}


results_ratio <- results_ratio %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))



# Plot coefficients
rdd_binary_coefs <- results_binary %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(0,0.15) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11)) +
  scale_y_continuous(labels=percent, breaks = seq(0,0.1,0.05), limits = c(0,0.13)) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )



ymin_binary <- 0
ymax_binary <- 0.75

rdd_binary_rdplot <- rdplot(
  y = post_data$relative_binary,
  x = post_data$margin,
  c = 0,
  p = 2,
  covs = cbind(post_data$party_name, post_data$year),
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  title = "",
  x.label = "Margin",
  y.label = "Percentage of respondents"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin_binary, ymax_binary), expand = c(0,0), labels=percent) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("Own voters more important\n(0/1)") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_binary$bw_est[1],
    xmax = results_binary$bw_est[1],
    ymin = ymin_binary,
    ymax = ymax_binary
  ),
  fill = "grey",
  alpha = 0.5) 



# Plot coefficients
rdd_ratio_coefs <- results_ratio %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11)) +
  scale_y_continuous(breaks = seq(0,0.1,0.05), limits = c(0,0.12)) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )



ymin_ratio <- -0.16
ymax_ratio <- 0.26


rdd_ratio_rdplot <- rdplot(
  y = post_data$relative_ratio,
  x = post_data$margin,
  c = 0,
  p = 2,
  covs = cbind(post_data$party_name, post_data$year),
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  title = "",
  x.label = "Margin",
  y.label = "(VOTERS-ALL)/max(VOTERS, ALL)"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin_ratio, ymax_ratio), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("Normalized relative\nimportance") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_binary$bw_est[1],
    xmax = results_binary$bw_est[1],
    ymin = ymin_ratio,
    ymax = ymax_ratio
  ),
  fill = "grey",
  alpha = 0.5) 


# Combine plots
plot_grid(
  rdd_binary_rdplot,
  rdd_ratio_rdplot,
  rdd_binary_coefs,
  rdd_ratio_coefs,
  
  align = "v",
  nrow = 2 ,   
  rel_heights = c(2.2/3, 1/3)
)



ggsave("plots/figure_d1.1.pdf", width = 1800, height = 1350, units = "px")






## Figure D2.1: Robustness check: LATEs with 95% confidence interval for different bandwidth specifications --------

# Set different bandwidths (ranging between 0.25 and 10, in 0.25 steps)
bandwidth <- seq(0.25, 10, 0.25)

# Empty frames to pipe in the results
results_own_40 <- results_all_40 <- results_relative_40 <- matrix(NA, ncol = 4, nrow = length(bandwidth))

# Change colnames of the empty data frames
colnames(results_own_40) <-
  colnames(results_all_40) <-
  colnames(results_relative_40) <-
  c("bandwidth",
    "LATE",
    "lower",
    "upper")


# Loop running the same estimation for the different bandwidths
for (i in 1:length(bandwidth)) {
  
  # RDD: Represent own voters
  rdd_own_bw <- rdrobust::rdrobust(y = post_data$rep_voters,
                                   x = post_data$margin,
                                   c = 0,
                                   p = 1,
                                   h = bandwidth[i],
                                   covs = cbind(post_data$party_name, post_data$year),
                                   vce = "nn",
                                   cluster = post_data$cluster,
                                   all = TRUE
  )
  
  # RDD: Represent all citizens in constituency
  rdd_all_bw <- rdrobust::rdrobust(y = post_data$rep_all,
                                   x = post_data$margin,
                                   c = 0,
                                   p = 1,
                                   h = bandwidth[i],
                                   covs = cbind(post_data$party_name, post_data$year),
                                   vce = "nn",
                                   cluster = post_data$cluster,
                                   all = TRUE
  )
  
  # RDD: Relative importance
  rdd_relative_bw <- rdrobust::rdrobust(y = post_data$relative_importance,
                                        x = post_data$margin,
                                        c = 0,
                                        p = 1,
                                        h = bandwidth[i],
                                        covs = cbind(post_data$party_name, post_data$year),
                                        vce = "nn",
                                        cluster = post_data$cluster,
                                        all = TRUE
  )
  
  
  # Piping in the results
  results_own_40[i,1] <- bandwidth[i]
  results_own_40[i,2] <- rdd_own_bw$Estimate[1]
  results_own_40[i,3] <- rdd_own_bw$ci[1,1]
  results_own_40[i,4] <- rdd_own_bw$ci[1,2]
  
  results_all_40[i,1] <- bandwidth[i]
  results_all_40[i,2] <- rdd_all_bw$Estimate[1]
  results_all_40[i,3] <- rdd_all_bw$ci[1,1]
  results_all_40[i,4] <- rdd_all_bw$ci[1,2]
  
  results_relative_40[i,1] <- bandwidth[i]
  results_relative_40[i,2] <- rdd_relative_bw$Estimate[1]
  results_relative_40[i,3] <- rdd_relative_bw$ci[1,1]
  results_relative_40[i,4] <- rdd_relative_bw$ci[1,2]
  
  
}


# Plot the results 

# Represent own voters
plot_40_own <- results_own_40 %>% 
  as.data.frame() %>% 
  ggplot(aes(x = bandwidth, y = LATE)) +
  theme_bw() +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.1) +
  geom_point(size = 0.7, position = position_dodge(width = 0.75)) +
  scale_y_continuous(limits = c(-0.96, 0.96), expand = c(0,0), breaks = seq(-0.75,0.75,0.25)) +
  scale_x_continuous(limits = c(-0.25,10.5), expand = c(0,0)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray20") +
  ggtitle("Should represent\nown voters\nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 12))+
  xlab(element_text("Bandwidth")) +
  theme(panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  )


# Represent all citizens in constituency
plot_40_all <- results_all_40 %>% 
  as.data.frame() %>% 
  ggplot(aes(x = bandwidth, y = LATE)) +
  theme_bw() + 
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.1) +
  geom_point(size = 0.7, position = position_dodge(width = 0.75)) +
  scale_y_continuous(limits = c(-0.96, 0.96), expand = c(0,0), breaks = seq(-0.75,0.75,0.25)) +
  scale_x_continuous(limits = c(-0.25,10.5), expand = c(0,0)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray20") +
  ggtitle("Should represent\nall citizens\nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 12))+
  xlab(element_text("Bandwidth")) +
  theme(panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  )


# Relative importance
plot_40_relative <- results_relative_40 %>% 
  as.data.frame() %>% 
  ggplot(aes(x = bandwidth, y = LATE)) +
  theme_bw() +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.1) +
  geom_point(size = 0.7, position = position_dodge(width = 0.75)) +
  scale_y_continuous(limits = c(-0.96, 0.96), expand = c(0,0), breaks = seq(-0.75,0.75,0.25)) +
  scale_x_continuous(limits = c(-0.25,10.5), expand = c(0,0)) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray20") +
  ggtitle("\nRelative\nimportance") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold", size = 12)) +
  xlab(element_text("Bandwidth")) +
  theme(panel.grid.minor.y = element_blank(),
        panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank()
  )


# Combine plots
plot_grid(
  plot_40_own,
  plot_40_all,
  plot_40_relative,
  
  align = "v",
  nrow = 1 
)


ggsave("plots/figure_d2.1.pdf", width = 1800, height = 900, units = "px")







## Figure D3.1: Robustness to different orders of the local-polynomial ---------------------------------------------

# Set range of polynomial orders
polynomials <- 1:4

specification <- c("Conventional", "Bias-corrected", "Robust")

# Empty frames to pipe in the results
results_own_polynomial <- results_all_polynomial <- results_relative_polynomial <- array(NA, dim = c(length(polynomials), 5, length(specification)))

colnames(results_own_polynomial) <-
  colnames(results_all_polynomial) <-
  colnames(results_relative_polynomial) <-
  c("polynomial",
    "LATE",
    "lower",
    "upper",
    "specification"
  )



# Loop running the same estimation for the different bandwidths
for (i in 1:length(polynomials)) {
  
  # RDD: Represent own voters
  rdd_own_polynomial <- rdrobust::rdrobust(y = post_data$rep_voters,
                                           x = post_data$margin,
                                           c = 0,
                                           p = polynomials[i],
                                           covs = cbind(post_data$party_name, post_data$year),
                                           vce = "nn",
                                           cluster = post_data$cluster,
                                           all = TRUE
  )
  
  # RDD: Represent all citizens in constituency
  rdd_all_polynomial <- rdrobust::rdrobust(y = post_data$rep_all,
                                           x = post_data$margin,
                                           c = 0,
                                           p = polynomials[i],
                                           covs = cbind(post_data$party_name, post_data$year),
                                           vce = "nn",
                                           cluster = post_data$cluster,
                                           all = TRUE
  )
  
  # RDD: Relative importance
  rdd_relative_polynomial <- rdrobust::rdrobust(y = post_data$relative_importance,
                                                x = post_data$margin,
                                                c = 0,
                                                p = polynomials[i],
                                                covs = cbind(post_data$party_name, post_data$year),
                                                vce = "nn",
                                                cluster = post_data$cluster,
                                                all = TRUE
  )
  
  
  
  
  
  for (j in seq_along(specification)) {
    
    
    temp <- dplyr::if_else(j == 1, "tau.us", "tau.bc")
    
    results_own_polynomial[i,1,j] <- polynomials[i]
    results_own_polynomial[i,2,j] <- rdd_own_polynomial$Estimate[,temp]
    results_own_polynomial[i,3,j] <- rdd_own_polynomial$ci[j,1]
    results_own_polynomial[i,4,j] <- rdd_own_polynomial$ci[j,2]
    results_own_polynomial[i,5,j] <- specification[j]
    
    
    results_all_polynomial[i,1,j] <- polynomials[i]
    results_all_polynomial[i,2,j] <- rdd_all_polynomial$Estimate[,temp]
    results_all_polynomial[i,3,j] <- rdd_all_polynomial$ci[j,1]
    results_all_polynomial[i,4,j] <- rdd_all_polynomial$ci[j,2]
    results_all_polynomial[i,5,j] <- specification[j]
    
    results_relative_polynomial[i,1,j] <- polynomials[i]
    results_relative_polynomial[i,2,j] <- rdd_relative_polynomial$Estimate[,temp]
    results_relative_polynomial[i,3,j] <- rdd_relative_polynomial$ci[j,1]
    results_relative_polynomial[i,4,j] <- rdd_relative_polynomial$ci[j,2]
    results_relative_polynomial[i,5,j] <- specification[j]
    
  }
  
}


temp_own <- rbind(results_own_polynomial[,,1],
                  results_own_polynomial[,,2],
                  results_own_polynomial[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-("specification")), as.numeric) %>% 
  mutate(dv = "Own voters in constituency")

temp_all <- rbind(results_all_polynomial[,,1],
                  results_all_polynomial[,,2],
                  results_all_polynomial[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-("specification")), as.numeric) %>% 
  mutate(dv = "All citizens in constituency")

temp_relative <- rbind(results_relative_polynomial[,,1],
                       results_relative_polynomial[,,2],
                       results_relative_polynomial[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-("specification")), as.numeric) %>% 
  mutate(dv = "Relative importance")


results_polynomial <- rbind(temp_own, temp_all, temp_relative) %>% 
  mutate(specification = factor(specification,
                                c("Robust", "Bias-corrected", "Conventional"),
                                labels = c("Robust", "Bias-corrected", "Conventional")
  )) %>%
  mutate(dv = factor(dv,
                     c("Own voters in constituency", "All citizens in constituency", "Relative importance"),
                     labels = c("Should represent own\nvoters in constituency", "Should represent all\ncitizens in constituency", "Relative\nimportance")
  )) %>% 
  mutate(polynomial = factor(polynomial,
                             c(1:4),
                             labels = c("1", "2", "3", "4")
  ))



# Plot results
results_polynomial %>% 
  ggplot(aes(LATE, polynomial)) +
  theme_bw() +
  geom_pointrange(
    aes(xmin = lower,
        xmax = upper,
        color = specification,
        shape = specification
    ),
    position = position_dodge(0.5),
    lwd = 0.6
  ) +
  scale_y_discrete(limits=rev) +
  facet_grid(cols = vars(dv), scales = "free_x") +
  ggh4x::facetted_pos_scales(
    x = list(
      dv == "Should represent own\nvoters in constituency" ~ scale_x_continuous(limits = c(-0.1,0.4), breaks = round(seq(-0.1,0.4,0.1), digits = 1)),
      dv == "Should represent all\ncitizens in constituency" ~ scale_x_continuous(limits = c(-0.5,0), breaks = round(seq(-0.5,0,0.1), digits = 1)),
      dv == "Relative\nimportance" ~ scale_x_continuous(limits = c(-0.05,0.65), breaks = seq(0,0.65, 0.2) )
    )
  ) +
  geom_vline(xintercept = 0,
             lty = 2,
             colour = gray(1 / 2)) +
  ylab(element_text("Polynomial (conventional)")) +
  theme(legend.title=element_blank()) +
  scale_fill_discrete(breaks=c("Conventional", "Bias-corrected", "Robust")) +
  scale_colour_manual(name="",
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  theme(legend.position="bottom") +
  # Reverse order of legend
  guides(color = guide_legend(reverse=TRUE),
         shape = guide_legend(reverse=TRUE)
  ) + 
  theme(panel.grid.minor.x = element_blank()) +
  
  theme(strip.text.x = element_text(size = 11, face = "bold"))


ggsave("plots/figure_d3.1.pdf", width = 1950, height = 1350, units = "px")




## Figure D4.1: Additional robustness checks -------------------------------

# Recalculate main findings 

rdd_own <- rdrobust::rdrobust(y = post_data$rep_voters, 
                              x = post_data$margin, 
                              c = 0, 
                              p = 1,
                              q = 2,
                              covs = cbind(post_data$party_name, post_data$year), 
                              vce = "nn", 
                              cluster = post_data$cluster, 
                              all = TRUE
)

rdd_all <- rdrobust::rdrobust(y = post_data$rep_all, 
                              x = post_data$margin, 
                              c = 0, 
                              p = 1,
                              q = 2,
                              covs = cbind(post_data$party_name, post_data$year), 
                              vce = "nn", 
                              cluster = post_data$cluster, 
                              all = TRUE
)

rdd_relative <- rdrobust::rdrobust(y = post_data$relative_importance, 
                                   x = post_data$margin, 
                                   c = 0, 
                                   p = 1,
                                   q = 2,
                                   covs = cbind(post_data$party_name, post_data$year), 
                                   vce = "nn", 
                                   cluster = post_data$cluster, 
                                   all = TRUE
)



# With socio-demographic weights 
rdd_own_weights <- rdrobust::rdrobust(y = post_data$rep_voters, 
                                      x = post_data$margin, 
                                      c = 0, 
                                      p = 1,
                                      covs = cbind(post_data$party_name, post_data$year), 
                                      vce = "nn", 
                                      cluster = post_data$cluster, 
                                      weights = post_data$weight,
                                      all = TRUE
)

rdd_all_weights <- rdrobust::rdrobust(y = post_data$rep_all, 
                                      x = post_data$margin, 
                                      c = 0, 
                                      p = 1,
                                      covs = cbind(post_data$party_name, post_data$year), 
                                      vce = "nn", 
                                      cluster = post_data$cluster, 
                                      weights = post_data$weight,
                                      all = TRUE
)

rdd_relative_weights <- rdrobust::rdrobust(y = post_data$relative_importance, 
                                           x = post_data$margin, 
                                           c = 0, 
                                           p = 1,
                                           covs = cbind(post_data$party_name, post_data$year), 
                                           vce = "nn", 
                                           cluster = post_data$cluster, 
                                           weights = post_data$weight,
                                           all = TRUE
)



# Only electoral races in which all six major parties ran (CDU/CSU, SPD, Greens, FDP, Left, AfD)
post_data_party6 <- post_data %>% filter(n_parties == 6)

rdd_own_party6 <- rdrobust::rdrobust(y = post_data_party6$rep_voters, 
                                     x = post_data_party6$margin, 
                                     c = 0, 
                                     p = 1,
                                     q = 2,
                                     vce = "nn", 
                                     covs = cbind(post_data_party6$party_name, post_data_party6$year),
                                     cluster = post_data_party6$cluster,
                                     all = TRUE
)

rdd_all_party6 <- rdrobust::rdrobust(y = post_data_party6$rep_all, 
                                     x = post_data_party6$margin, 
                                     c = 0, 
                                     p = 1,
                                     q = 2,
                                     vce = "nn", 
                                     covs = cbind(post_data_party6$party_name, post_data_party6$year),
                                     cluster = post_data_party6$cluster, 
                                     all = TRUE
)

rdd_relative_party6 <- rdrobust::rdrobust(y = post_data_party6$relative_importance, 
                                          x = post_data_party6$margin, 
                                          c = 0, 
                                          p = 1,
                                          q = 2,
                                          vce = "nn", 
                                          covs = cbind(post_data_party6$party_name, post_data_party6$year),
                                          cluster = post_data_party6$cluster, 
                                          all = TRUE
)


# Excluding party + year dummies from the analyses
rdd_own_nodummy <- rdrobust::rdrobust(y = post_data$rep_voters, 
                                      x = post_data$margin, 
                                      c = 0, 
                                      p = 1,
                                      q = 2,
                                      vce = "nn", 
                                      cluster = post_data$cluster, 
                                      all = TRUE
)

rdd_all_nodummy <- rdrobust::rdrobust(y = post_data$rep_all, 
                                      x = post_data$margin, 
                                      c = 0, 
                                      p = 1,
                                      q = 2,
                                      vce = "nn",
                                      cluster = post_data$cluster, 
                                      all = TRUE
)

rdd_relative_nodummy <- rdrobust::rdrobust(y = post_data$relative_importance, 
                                           x = post_data$margin, 
                                           c = 0, 
                                           p = 1,
                                           q = 2,
                                           vce = "nn", 
                                           cluster = post_data$cluster, 
                                           all = TRUE
)










model_specification <- c("Conventional", "Bias-corrected", "Robust")

models <- c("rdd_own",
            "rdd_all",
            "rdd_relative",
            "rdd_own_weights",
            "rdd_all_weights",
            "rdd_relative_weights",
            "rdd_own_party6",
            "rdd_all_party6",
            "rdd_relative_party6",
            "rdd_own_nodummy",
            "rdd_all_nodummy",
            "rdd_relative_nodummy"
)

results_robustness <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_,
  model = NA_character_
)


results_robustness_final <- results_robustness

# Pipe in the results
for (m in seq_along(models)) {
  
  
  for(i in 1: length(model_specification)){
    
    rdd <- get(models[m])
    
    results_robustness[i, "specification"] <- model_specification[i]
    results_robustness[i, "LATE"] <- rdd$coef[i] 
    results_robustness[i, "lower"] <- rdd$ci[i,1]
    results_robustness[i, "upper"] <- rdd$ci[i,2]
    results_robustness[i, "SE"] <- rdd$se[i]
    results_robustness[i, "p"] <- rdd$pv[i]
    results_robustness[i, "bw_est"] <- rdd$bws[1,1]
    results_robustness[i, "bw_bias"] <- rdd$bws[2,1]
    results_robustness[i, "total_obs"] <- rdd$N %>% sum()
    results_robustness[i, "eff_obs"] <- if_else(i == 1, rdd$N_h %>% sum(), rdd$N_b %>% sum())
    results_robustness[i, "left"] <- if_else(i == 1, rdd$N_h[1], rdd$N_b[1])
    results_robustness[i, "right"] <- if_else(i == 1, rdd$N_h[2], rdd$N_b[2])
    results_robustness[i, "model"] <- models[m]
  }
  
  
  results_robustness_final <- rbind(results_robustness_final, results_robustness)
  
  
  
  
  
}  


# PLot the results
results_robustness_final %>% 
  filter(!is.na(specification)) %>% 
  mutate(robustness = case_when(
    grepl("weights", model) ~ "weights",
    grepl("party6", model) ~ "party6",
    grepl("nodummy", model) ~ "nodummy",
    TRUE ~ "main"
  )) %>% 
  
  mutate(robustness = factor(robustness,
                             levels = c("main",
                                        "weights",
                                        "party6",
                                        "nodummy"
                             ),
                             labels = c("Main\nfindings",
                                        "With socio-demo-\ngraphic weights",
                                        "Only races with all\nsix major parties",
                                        "Excluding party\n+ year dummies"
                             )
  )) %>% 
  
  mutate(dv = case_when(
    grepl("own", model) ~ "own",
    grepl("all", model) ~ "all",
    grepl("relative", model) ~ "relative"
  )) %>%  
  
  mutate(dv = factor(dv, 
                     levels = c("own", "all", "relative"),
                     labels = c("Should represent own\nvoters in constituency",
                                "Should represent all\ncitizens in constituency",
                                "Relative\nimportance"
                     )
  )) %>% 
  
  mutate(specification = factor(specification,
                                levels = c("Robust", "Bias-corrected", "Conventional")
  )) %>% 
  
  
  
  
  ggplot(aes(LATE, robustness)) +
  theme_bw() +
  geom_pointrange(
    aes(xmin = lower,
        xmax = upper,
        color = specification,
        shape = specification
    ),
    position = position_dodge(0.5),
    lwd = 0.6
  ) +
  scale_y_discrete(limits=rev) +
  facet_grid(cols = vars(dv), scales = "free_x") +
  geom_vline(xintercept = 0,
             lty = 2,
             colour = gray(1 / 2)) +
  ylab(element_blank()) +
  theme(legend.title=element_blank()) +
  scale_fill_discrete(breaks=c("Conventional", "Bias-corrected", "Robust")) +
  scale_colour_manual(name="",
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  theme(legend.position="bottom") +
  # Reverse order of legend
  guides(color = guide_legend(reverse=TRUE),
         shape = guide_legend(reverse=TRUE)
  ) + 
  theme(panel.grid.minor.x = element_blank()) +
  
  theme(strip.text.x = element_text(size = 11, face = "bold"))


ggsave("plots/figure_d4.1.pdf", width = 2150, height = 1350, units = "px")







## Figure E1.1: Balance statistics -----------------------------------------

# Create standardized age variable
post_data <- post_data %>%
  rowwise() %>%
  mutate(age_z = (age - mean(post_data$age, na.rm = TRUE)) / sd(post_data$age, na.rm = TRUE))

variables <- c("age_z", "male", "education", "interest", "knowledge")
model_specification <- c("Conventional", "Bias corrected", "Robust")



results_balance <- data.frame(
  dv = NA_character_,
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_
)

results_balance_post <- NULL



# Repeat RDD for different socio-demographic factors and political variables to check for balance at the cutoff
for (i in 1:length(variables)) {
  rdd <- rdrobust::rdrobust(y = post_data %>% pull(variables[i]),
                            x = post_data$margin,
                            c = 0,
                            p = 1,
                            q = 2,
                            vce = "nn",
                            cluster = cbind(post_data$cluster),
                            all = TRUE
  )
  
  for (j in 1:length(model_specification)) {
    
    results_balance[j, "dv"] <- variables[i]  
    results_balance[j, "specification"] <- model_specification[j]
    results_balance[j, "LATE"] <- rdd$coef[j] 
    results_balance[j, "lower"] <- rdd$ci[j,1]
    results_balance[j, "upper"] <- rdd$ci[j,2]
    
  }
  
  results_balance_post <- rbind(results_balance_post, results_balance)
  
  
}


# Plot the results
results_balance_post %>% 
  mutate(dv = fct_rev(factor(
    dv,
    levels = c("age_z", "male", "education", "interest", "knowledge"),
    labels = c("Age (z-standardized)", "Male (0/1)", "Education (1-5)", "Political interest (1-5)", "Political knowledge (0/1)")
  ))) %>% 
  mutate(specification = factor(
    specification,
    levels = c("Robust", "Bias corrected", "Conventional"),
    labels = c("Robust", "Bias corrected", "Conventional")
  )) %>%
  ggplot(aes(LATE, dv)) +
  theme_bw() +
  geom_pointrange(
    aes(xmin = lower,
        xmax = upper,
        color = specification,
        shape = specification
    ),
    position = position_dodge(0.7),
    lwd = 0.6
  ) +
  geom_vline(xintercept = 0,
             lty = 2,
             colour = gray(1 / 2)) +
  ylab(element_blank()) +
  theme(legend.title=element_blank()) +
  scale_fill_discrete(breaks=c('Conventional', 'Bias corrected', 'Robust')) +
  scale_colour_manual(name="",  
                      values = c("Conventional" = viridis(3)[1], "Bias corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias corrected" = 2, "Robust" = 15)) +
  # Reverse order of legend
  guides(color = guide_legend(reverse=TRUE),
         shape = guide_legend(reverse=TRUE)
  )

ggsave("plots/figure_e1.1.pdf", width = 1800, height = 900, units = "px")






## Figure E2.1: Manipulation checks for the forcing variable ---------------

## Run McCrary (2008) version of the test

# Get p-value of the test
p_val_mccrary <- DCdensity(runvar = post_data$margin, 
                           c = 0,
                           plot = FALSE
) %>% round(digits = 3)

# Get plot from the package
p1 <- ~{
  par(
    mar = c(2.5,2,1,0.5),
    mgp = c(2, 0.1, 0),
    cex.axis = 0.74,
    tck=-0.01,
    las=1
  )
  DCdensity(runvar = post_data$margin,
            c = 0,
            plot = TRUE
  )
  abline(v = 0, lty = 2) # Add vertical, dotted line at 0 (the cutoff in the RDD)
  text(x = 22, y = 0.028, paste0("p = ", p_val_mccrary)) # Add p-value to the plot
}



# Cattaneo et al. (2020) maniupaltion check

# Run the model
rdd <- rddensity(X = post_data$margin, vce = "jackknife")

# Get the plot
cattaneo_plot <- rddensity::rdplotdensity(rdd,
                                          X = post_data$margin,
                                          plotRange = c(-40, 40)
) 

# Get the p-value
p_val_cattaneo <- rdd$bino$pval[5] %>% round(digits = 3)

# Finsih the plot and add p-value to it
p2 <- cattaneo_plot$Estplot +
  geom_vline(xintercept = 0, lty = 2) +
  scale_x_continuous(limits = c(-40,40), expand = c(0,0)) +
  scale_y_continuous(limits = c(0,0.0315), expand = c(0.03,0.0005), breaks = seq(0,0.03,0.005)) + 
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
        panel.background = element_blank(), axis.line = element_line(colour = "black")) +
  theme(axis.text.x=element_text(colour="black"),
        axis.text.y=element_text(colour="black")) +
  ggtitle(element_blank()) +
  annotate("text", x=22, y=0.028, label= paste0("p = ", p_val_cattaneo), size = 4.3)


# Combine the two plots together with titles
p <- plot_grid(p1, p2, labels = c('McCrary (2008)', 'Cattaneo et al. (2020)'), label_size = 10)

# Common x-axis label
x.grob <- textGrob("Margin of victory/loss", 
                   gp=gpar(fontsize=15))

# Common y-axis label
y.grob <- textGrob("Density", 
                   gp=gpar(fontsize=15), rot=90)

# Arrange plots
grid.arrange(arrangeGrob(p, left = y.grob, bottom = x.grob))

ggsave("plots/figure_e2.1.pdf", width = 1950, height = 1350, units = "px")









## Figure F1.1: Placebo test: LATEs for respondents interviewed before the elections ---------------------------------------

# NOTE: The following code follows the exact same structure as the code for Figure 2 (see above), just for respondents interviewed before the election dates (as placebo checks).

model_specification <- c("Conventional", "Bias corrected", "Robust")

# Represent own voters 
rdd_own_pre <- rdrobust::rdrobust(y = pre_data$rep_voters, 
                                  x = pre_data$margin, 
                                  c = 0, 
                                  p = 1,
                                  q = 2,
                                  covs = cbind(pre_data$party_name, pre_data$year), 
                                  vce = "nn", 
                                  cluster = pre_data$cluster, 
                                  all = TRUE
)

results_own_pre <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_own_pre[i, "specification"] <- model_specification[i]
  results_own_pre[i, "LATE"] <- rdd_own_pre$coef[i] 
  results_own_pre[i, "lower"] <- rdd_own_pre$ci[i,1]
  results_own_pre[i, "upper"] <- rdd_own_pre$ci[i,2]
  results_own_pre[i, "SE"] <- rdd_own_pre$se[i]
  results_own_pre[i, "p"] <- rdd_own_pre$pv[i]
  results_own_pre[i, "bw_est"] <- rdd_own_pre$bws[1,1]
  results_own_pre[i, "bw_bias"] <- rdd_own_pre$bws[2,1]
  results_own_pre[i, "total_obs"] <- rdd_own_pre$N %>% sum()
  results_own_pre[i, "eff_obs"] <- if_else(i == 1, rdd_own_pre$N_h %>% sum(), rdd_own_pre$N_b %>% sum())
  results_own_pre[i, "left"] <- if_else(i == 1, rdd_own_pre$N_h[1], rdd_own_pre$N_b[1])
  results_own_pre[i, "right"] <- if_else(i == 1, rdd_own_pre$N_h[2], rdd_own_pre$N_b[2])
}



results_own_pre <- results_own_pre %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))



# Represent all citizens 
rdd_all_pre <- rdrobust::rdrobust(y = pre_data$rep_all, 
                                  x = pre_data$margin, 
                                  c = 0, 
                                  p = 1,
                                  q = 2,
                                  covs = cbind(pre_data$party_name, pre_data$year), 
                                  vce = "nn", 
                                  cluster = pre_data$cluster, 
                                  all = TRUE
)


model_specification <- c("Conventional", "Bias corrected", "Robust")

results_all_pre <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_all_pre[i, "specification"] <- model_specification[i]
  results_all_pre[i, "LATE"] <- rdd_all_pre$coef[i] 
  results_all_pre[i, "lower"] <- rdd_all_pre$ci[i,1]
  results_all_pre[i, "upper"] <- rdd_all_pre$ci[i,2]
  results_all_pre[i, "SE"] <- rdd_all_pre$se[i]
  results_all_pre[i, "p"] <- rdd_all_pre$pv[i]
  results_all_pre[i, "bw_est"] <- rdd_all_pre$bws[1,1]
  results_all_pre[i, "bw_bias"] <- rdd_all_pre$bws[2,1]
  results_all_pre[i, "total_obs"] <- rdd_all_pre$N %>% sum()
  results_all_pre[i, "eff_obs"] <- if_else(i == 1, rdd_all_pre$N_h %>% sum(), rdd_all_pre$N_b %>% sum())
  results_all_pre[i, "left"] <- if_else(i == 1, rdd_all_pre$N_h[1], rdd_all_pre$N_b[1])
  results_all_pre[i, "right"] <- if_else(i == 1, rdd_all_pre$N_h[2], rdd_all_pre$N_b[2])
}


results_all_pre <- results_all_pre %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))



# Relative importance
rdd_relative_pre <- rdrobust::rdrobust(y = pre_data$relative_importance, 
                                       x = pre_data$margin, 
                                       c = 0, 
                                       p = 1,
                                       q = 2,
                                       covs = cbind(pre_data$party_name, pre_data$year), 
                                       vce = "nn", 
                                       cluster = pre_data$cluster, 
                                       all = TRUE
)

model_specification <- c("Conventional", "Bias corrected", "Robust")

results_relative_pre <- data.frame(
  specification = NA_character_,
  LATE = NA_real_,
  lower = NA_real_,
  upper = NA_real_,
  SE = NA_real_,
  p = NA_real_,
  bw_est = NA_real_,
  bw_bias = NA_real_,
  total_obs = NA_real_,
  eff_obs = NA_real_,
  left = NA_real_,
  right = NA_real_
)


for(i in 1: length(model_specification)){
  
  results_relative_pre[i, "specification"] <- model_specification[i]
  results_relative_pre[i, "LATE"] <- rdd_relative_pre$coef[i] 
  results_relative_pre[i, "lower"] <- rdd_relative_pre$ci[i,1]
  results_relative_pre[i, "upper"] <- rdd_relative_pre$ci[i,2]
  results_relative_pre[i, "SE"] <- rdd_relative_pre$se[i]
  results_relative_pre[i, "p"] <- rdd_relative_pre$pv[i]
  results_relative_pre[i, "bw_est"] <- rdd_relative_pre$bws[1,1]
  results_relative_pre[i, "bw_bias"] <- rdd_relative_pre$bws[2,1]
  results_relative_pre[i, "total_obs"] <- rdd_relative_pre$N %>% sum()
  results_relative_pre[i, "eff_obs"] <- if_else(i == 1, rdd_relative_pre$N_h %>% sum(), rdd_relative_pre$N_b %>% sum())
  results_relative_pre[i, "left"] <- if_else(i == 1, rdd_relative_pre$N_h[1], rdd_relative_pre$N_b[1])
  results_relative_pre[i, "right"] <- if_else(i == 1, rdd_relative_pre$N_h[2], rdd_relative_pre$N_b[2])
}



results_relative_pre <- results_relative_pre %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias corrected", "Robust"),
                                labels = c("Conven-\ntional", "Bias \ncorrec.", "Robust")
  ))






### Own voters
# Plot the coefficients
rdd_own_coefs_pre <- results_own_pre %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(-0.3,0.5) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11)) +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )
# RDD plot:

ymin <- 3
ymax <- 5

rdd_own_rdplot_pre <- rdplot(
  y = pre_data$rep_voters,
  x = pre_data$margin,
  c = 0,
  covs = cbind(pre_data$party_name, pre_data$year),
  #weights = pre_data$weight,
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  p = 2,
  title = "",
  x.label = "Margin",
  y.label = "Importance"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin, ymax), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("Should represent\nown voters \nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_own_pre$bw_est[1],
    xmax = results_own_pre$bw_est[1],
    ymin = ymin,
    ymax = ymax
  ),
  fill = "grey",
  alpha = 0.5)

### All voters
# Plot the coefficients
rdd_all_coefs_pre <- results_all_pre %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(-0.3,0.5) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11))+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )
# RDD plot:

ymin <- 3
ymax <- 5

rdd_all_rdplot_pre <- rdplot(
  y = pre_data$rep_all,
  x = pre_data$margin,
  c = 0,
  covs = cbind(pre_data$party_name, pre_data$year),
  #weights = pre_data$weight,
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  p = 2,
  title = "",
  x.label = "Margin",
  y.label = "Importance"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin, ymax), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("Should represent\nall citizens \nin constituency") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_all_pre$bw_est[1],
    xmax = results_all_pre$bw_est[1],
    ymin = ymin,
    ymax = ymax
  ),
  fill = "grey",
  alpha = 0.5)

### Relative importance
# Plot the coefficients
rdd_relative_coefs_pre <- results_relative_pre %>%
  ggplot(aes(specification, LATE)) + 
  geom_hline(aes(yintercept = 0), linetype = 2, color = "gray20") + 
  ylim(-0.3,0.5) +
  geom_linerange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.75), size = 0.5) + 
  geom_point(size = 1.75, position = position_dodge(width = 0.75)) + 
  theme_minimal( ) +
  xlab(element_blank()) + 
  ylab("LATE") +
  theme(axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 11),
        axis.title = element_text(size = 11))+
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.y = element_blank()
  )
# RDD plot:

ymin_relative <- -0.5
ymax_relative <- 1

rdd_relative_rdplot_pre <- rdplot(
  y = pre_data$relative_importance,
  x = pre_data$margin,
  c = 0,
  covs = cbind(pre_data$party_name, pre_data$year),
  #weights = pre_data$weight,
  x.lim = c(-29, 29),
  #nbins = c(20,20),
  #ci = 95,
  #shade = TRUE,
  p = 2,
  title = "",
  x.label = "Margin",
  y.label = "Relative importance"
)$rdplot +
  scale_x_continuous(breaks = seq(-30, 30, 10)) +
  scale_y_continuous(limits = c(ymin_relative, ymax_relative), expand = c(0,0)) +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.background = element_blank(),
    axis.title = element_text(size = 11),
    axis.text.x = element_text(size = 11),
    axis.text.y = element_text(size = 11)
  ) +
  ggtitle("\nRelative \nimportance") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  # Add visualization of bandwidths
  geom_rect(aes(
    xmin = -results_relative_pre$bw_est[1],
    xmax = results_relative_pre$bw_est[1],
    ymin = ymin_relative,
    ymax = ymax_relative
  ),
  fill = "grey",
  alpha = 0.5)


# Combine all plots

plot_grid(
  rdd_own_rdplot_pre,
  rdd_all_rdplot_pre,
  rdd_relative_rdplot_pre,
  rdd_own_coefs_pre,
  rdd_all_coefs_pre,
  rdd_relative_coefs_pre,
  
  align = "v",
  nrow = 2 ,   
  rel_heights = c(2.3/3, 1/3)
)



ggsave("plots/figure_f1.1.pdf", width = 2100, height = 1250, units = "px")






## Figure F2.1: LATEs for several placebo cutoffs --------------------------

# NOTE: The code follows the exacts same structure as for Figure D3.1, just that here the cutoff points are varied.

placebo_cutoff <- seq(-15, 15, 2.5)

results_placebo_conventional <- results_placebo_bias <- results_placebo_robust <- matrix(NA, ncol = 5, nrow = length(placebo_cutoff))

colnames(results_placebo_conventional) <- colnames(results_placebo_bias) <- colnames(results_placebo_robust) <-
  c("placebo",
    "LATE",
    "lower",
    "upper",
    "specification"
  )


for (i in 1:length(placebo_cutoff)) {
  
  rdd_relative_placebo <- rdrobust::rdrobust(y = post_data$relative_importance,
                                             x = post_data$margin,
                                             c = placebo_cutoff[i],
                                             p = 1,
                                             covs = cbind(post_data$party_name, post_data$year),
                                             cluster = post_data$cluster,
                                             vce = "nn",
                                             all = TRUE
  )
  
  results_placebo_conventional[i,1] <- placebo_cutoff[i]
  results_placebo_conventional[i,2] <- rdd_relative_placebo$Estimate[1]
  results_placebo_conventional[i,3] <- rdd_relative_placebo$ci[1,1]
  results_placebo_conventional[i,4] <- rdd_relative_placebo$ci[1,2]
  results_placebo_conventional[i,5] <- "Conventional"
  
  results_placebo_bias[i,1] <- placebo_cutoff[i]
  results_placebo_bias[i,2] <- rdd_relative_placebo$Estimate[2]
  results_placebo_bias[i,3] <- rdd_relative_placebo$ci[2,1]
  results_placebo_bias[i,4] <- rdd_relative_placebo$ci[2,2]
  results_placebo_bias[i,5] <- "Bias-corrected"
  
  results_placebo_robust[i,1] <- placebo_cutoff[i]
  results_placebo_robust[i,2] <- rdd_relative_placebo$Estimate[2]
  results_placebo_robust[i,3] <- rdd_relative_placebo$ci[3,1]
  results_placebo_robust[i,4] <- rdd_relative_placebo$ci[3,2]
  results_placebo_robust[i,5] <- "Robust"
  
  
}

results_relative_placebo <- rbind(results_placebo_conventional, 
                                  results_placebo_bias,
                                  results_placebo_robust
) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-c("specification", "placebo")), as.numeric) %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias-corrected", "Robust"),
                                labels = c("Conventional", "Bias-corrected", "Robust")
  )) %>% 
  mutate(placebo = factor(placebo,
                          as.character(placebo_cutoff),
                          labels = as.character(placebo_cutoff)
  ))





# With sociodemographic weights

results_placebo_conventional_weights <- results_placebo_bias_weights <- results_placebo_robust_weights <- matrix(NA, ncol = 5, nrow = length(placebo_cutoff))

colnames(results_placebo_conventional_weights) <- colnames(results_placebo_bias_weights) <- colnames(results_placebo_robust_weights) <-
  c("placebo",
    "LATE",
    "lower",
    "upper",
    "specification"
  )


for (i in 1:length(placebo_cutoff)) {
  
  rdd_relative_placebo <- rdrobust::rdrobust(y = post_data$relative_importance,
                                             x = post_data$margin,
                                             c = placebo_cutoff[i],
                                             p = 1, 
                                             covs = cbind(post_data$party_name, post_data$year),
                                             cluster = post_data$cluster,
                                             weights = post_data$weight,
                                             vce = "nn",
                                             all = TRUE
  )
  
  results_placebo_conventional_weights[i,1] <- placebo_cutoff[i]
  results_placebo_conventional_weights[i,2] <- rdd_relative_placebo$Estimate[1]
  results_placebo_conventional_weights[i,3] <- rdd_relative_placebo$ci[1,1]
  results_placebo_conventional_weights[i,4] <- rdd_relative_placebo$ci[1,2]
  results_placebo_conventional_weights[i,5] <- "Conventional"
  
  results_placebo_bias_weights[i,1] <- placebo_cutoff[i]
  results_placebo_bias_weights[i,2] <- rdd_relative_placebo$Estimate[2]
  results_placebo_bias_weights[i,3] <- rdd_relative_placebo$ci[2,1]
  results_placebo_bias_weights[i,4] <- rdd_relative_placebo$ci[2,2]
  results_placebo_bias_weights[i,5] <- "Bias-corrected"
  
  results_placebo_robust_weights[i,1] <- placebo_cutoff[i]
  results_placebo_robust_weights[i,2] <- rdd_relative_placebo$Estimate[2]
  results_placebo_robust_weights[i,3] <- rdd_relative_placebo$ci[3,1]
  results_placebo_robust_weights[i,4] <- rdd_relative_placebo$ci[3,2]
  results_placebo_robust_weights[i,5] <- "Robust"
  
  
}

results_relative_placebo_weights <- rbind(results_placebo_conventional_weights, 
                                          results_placebo_bias_weights,
                                          results_placebo_robust_weights
) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-c("specification", "placebo")), as.numeric) %>% 
  mutate(specification = factor(specification,
                                c("Conventional", "Bias-corrected", "Robust"),
                                labels = c("Conventional", "Bias-corrected", "Robust")
  )) %>% 
  mutate(placebo = factor(placebo,
                          as.character(placebo_cutoff),
                          labels = as.character(placebo_cutoff)
  ))


plot_placebo_relative <- results_relative_placebo %>%
  as.data.frame() %>%
  ggplot(aes(x = placebo, y = LATE, col = specification, shape = specification)) +
  geom_pointrange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.5)) +
  ylim(-0.6, 0.6) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray20") +
  theme_minimal( ) +
  xlab(element_text("Placebo cutoff")) +
  ggtitle("Placebo tests: without sociodemographic weights") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  scale_fill_discrete(breaks=c('Conventional', 'Bias-corrected', 'Robust')) +
  scale_colour_manual(name="",  
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  
  
  theme(legend.position="none")

placebo_plot <- plot_placebo_relative + 
  annotate(
    "rect",
    xmin = 6.5,
    xmax = 7.5,
    ymin = -0.6,
    ymax = 0.6,
    alpha = 0.1,
    fill = "red"
  ) +
  annotate(
    "text",
    x = 7,
    y = -0.45,
    label = "Actual \ncutoff",
    size = 3.5
  )





plot_placebo_relative_weights <- results_relative_placebo_weights %>%
  as.data.frame() %>%
  ggplot(aes(x = placebo, y = LATE, col = specification, shape = specification)) +
  geom_pointrange(aes(ymin = lower, ymax = upper), position = position_dodge(width = 0.5)) +
  ylim(-0.6, 0.6) +
  geom_hline(yintercept = 0, linetype = 2, color = "gray20") +
  theme_minimal( ) +
  xlab(element_text("Placebo cutoff")) +
  ggtitle("Placebo tests: with sociodemographic weights") +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  scale_fill_discrete(breaks=c('Conventional', 'Bias-corrected', 'Robust')) +
  scale_colour_manual(name="",  
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  
  
  theme(legend.position="bottom")

placebo_plot_weights <- plot_placebo_relative_weights + 
  annotate(
    "rect",
    xmin = 6.5,
    xmax = 7.5,
    ymin = -0.6,
    ymax = 0.6,
    alpha = 0.1,
    fill = "red"
  ) +
  annotate(
    "text",
    x = 7,
    y = -0.45,
    label = "Actual \ncutoff",
    size = 3.5
  )



legend <- get_plot_component(placebo_plot_weights, 'guide-box-bottom', return_all = TRUE)


cowplot::plot_grid(placebo_plot + theme(legend.position="none"),
                   placebo_plot_weights + theme(legend.position="none"),
                   legend,
                   nrow = 3,
                   rel_heights = c(1,1,0.2)
)

ggsave("plots/figure_f2.1.pdf", width = 1800, height = 1800, units = "px")





## Figure G1.1: Changes in mean perceived relative importance between respondents interviewed before and after the election --------

# Use data set that includes both respondents interviewed before and those after the elections
representation_data_reduced <- representation_data %>% 
  
  # Remove missing values 
  filter(!is.na(relative_importance)) %>% 
  
  # Only take respondents whose preferred candidate won or lost by a margin smaller than 2 percentage points
  filter(abs(margin) <= 2) %>% 
  
  # Split sample into winner and losers
  mutate(vote_winner = dplyr::if_else(margin > 0, 1, 0)) %>% 
  
  # For each group: get mean values as well as 95% confidence intervals
  group_by(vote_winner, post) %>% 
  summarise(relative = mean(relative_importance),
            relative_sd = sd(relative_importance),
            n = n(),
            se = relative_sd/sqrt(n),
            lower = relative + qt(p = 0.025, df = n - 1) * se,
            upper = relative + qt(p = 0.975, df = n - 1) * se
  ) %>% 
  mutate(winner = dplyr::if_else(vote_winner == 1, "Winners", "Losers"),
         post = dplyr::if_else(post == 1, "Post-election", "Pre-election"),
         post = factor(post, 
                       levels = c("Pre-election", "Post-election")
         )
  )

# PLot the data
representation_data_reduced %>% 
  ggplot(aes(x = post, y = relative, color = winner, shape = winner)) +
  theme_bw() +
  geom_pointrange(aes(ymin = lower, ymax = upper), size = 1) + 
  geom_line(aes(group = winner)) +
  xlab(element_blank()) +
  ylab("Relative importance (mean)") +
  theme(legend.title = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(size = 12)
  ) +
  theme(axis.text.x = element_text(size = 12),
        axis.text.y = element_text(size = 12)) +
  scale_x_discrete(expand = c(0,0.25)) +
  scale_fill_discrete(breaks=c('Losers', 'Winners')) +
  scale_colour_manual(name="",  
                      values = c("Losers" = viridis(3)[2], "Winners" = viridis(3)[1])) +
  
  scale_shape_manual(name="",  
                     values = c("Losers" = 16, "Winners" = 15)) 


ggsave("plots/figure_g1.1.pdf", width = 1400, height = 1200, units = "px")






## Figure G2.1: RDD estimates by party affiliation of preferred can --------

# Specify the four relevant parties
parties <- c("cdu", "spd", "green", "afd")

# Set names of different model specifications 
specification <- c("Conventional", "Bias-corrected", "Robust")

# Empty arrays to pipe in the results
results_own_party <- results_all_party <- results_relative_party <- array(NA, dim = c(length(parties), 5, length(specification)))

colnames(results_own_party) <-
  colnames(results_all_party) <-
  colnames(results_relative_party) <-
  c("party",
    "LATE",
    "lower",
    "upper",
    "specification"
  )



# Loop for each party affiliation and outcome measures
for (i in 1:length(parties)) {
  
  # Only select respondents of the respective party
  data_party <- post_data %>% filter(party_name == parties[i])
  
  # RDD: Represent own voters
  rdd_own_party <- rdrobust::rdrobust(y = data_party$rep_voters,
                                      x = data_party$margin,
                                      c = 0,
                                      p = 1,
                                      q = 2,
                                      vce = "nn",
                                      cluster = data_party$cluster,
                                      covs = cbind(data_party$year),
                                      all = TRUE
  )
  
  # RDD: Represent all citizens in constituency 
  rdd_all_party <- rdrobust::rdrobust(y = data_party$rep_all,
                                      x = data_party$margin,
                                      c = 0,
                                      p = 1,
                                      q = 2,
                                      vce = "nn",
                                      cluster = data_party$cluster,
                                      covs = cbind(data_party$year),
                                      all = TRUE
  )
  
  # RDD: Relative importance
  rdd_relative_party <- rdrobust::rdrobust(y = data_party$relative_importance,
                                           x = data_party$margin,
                                           c = 0,
                                           p = 1,
                                           q = 2,
                                           vce = "nn",
                                           cluster = data_party$cluster,
                                           covs = cbind(data_party$year),
                                           all = TRUE
  )
  
  
  
  
  
  for (j in seq_along(specification)) {
    
    
    temp <- dplyr::if_else(j == 1, "tau.us", "tau.bc")
    
    results_own_party[i,1,j] <- parties[i]
    results_own_party[i,2,j] <- rdd_own_party$Estimate[,temp]
    results_own_party[i,3,j] <- rdd_own_party$ci[j,1]
    results_own_party[i,4,j] <- rdd_own_party$ci[j,2]
    results_own_party[i,5,j] <- specification[j]
    
    
    results_all_party[i,1,j] <- parties[i]
    results_all_party[i,2,j] <- rdd_all_party$Estimate[,temp]
    results_all_party[i,3,j] <- rdd_all_party$ci[j,1]
    results_all_party[i,4,j] <- rdd_all_party$ci[j,2]
    results_all_party[i,5,j] <- specification[j]
    
    results_relative_party[i,1,j] <- parties[i]
    results_relative_party[i,2,j] <- rdd_relative_party$Estimate[,temp]
    results_relative_party[i,3,j] <- rdd_relative_party$ci[j,1]
    results_relative_party[i,4,j] <- rdd_relative_party$ci[j,2]
    results_relative_party[i,5,j] <- specification[j]
    
  }
  
}





# Post-process these results 
temp_own <- rbind(results_own_party[,,1],
                  results_own_party[,,2],
                  results_own_party[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-c("specification", "party")), as.numeric) %>% 
  mutate(dv = "Own voters in constituency")

temp_all <- rbind(results_all_party[,,1],
                  results_all_party[,,2],
                  results_all_party[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-c("specification", "party")), as.numeric) %>% 
  mutate(dv = "All citizens in constituency")

temp_relative <- rbind(results_relative_party[,,1],
                       results_relative_party[,,2],
                       results_relative_party[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-c("specification", "party")), as.numeric) %>% 
  mutate(dv = "Relative importance")

# Bind data together
results_party <- rbind(temp_own, temp_all, temp_relative) %>% 
  mutate(specification = factor(specification,
                                c("Robust", "Bias-corrected", "Conventional"),
                                labels = c("Robust", "Bias-corrected", "Conventional")
  )) %>%
  mutate(dv = factor(dv,
                     c("Own voters in constituency", "All citizens in constituency", "Relative importance"),
                     labels = c("Should represent own\nvoters in constituency", "Should represent all\ncitizens in constituency", "Relative\nimportance")
  )) %>% 
  mutate(party = factor(party, 
                        levels = c("cdu", "spd", "green", "afd")
  ))


### Panel A --------------------------------------------

# Only results for CDU/CSU and SPD
results_party %>% 
  filter(party %in% c("cdu", "spd") ) %>% 
  mutate(party = factor(party, levels = c("cdu", "spd"), labels = c("CDU/CSU", "SPD"))) %>% 
  
  ggplot(aes(LATE, party)) +
  theme_bw() +
  geom_pointrange(
    aes(xmin = lower,
        xmax = upper,
        color = specification,
        shape = specification
    ),
    position = position_dodge(0.5),
    lwd = 0.6
  ) +
  scale_y_discrete(limits=rev) +
  facet_grid(cols = vars(dv), scales = "free_x") +
  ggh4x::facetted_pos_scales(
    x = list(
      dv == "Should represent own\nvoters in constituency" ~ scale_x_continuous(limits = c(-0.55,0.55), breaks = round(seq(-0.5,0.5,0.25), digits = 2)),
      dv == "Should represent all\ncitizens in constituency" ~ scale_x_continuous(limits = c(-0.55,0.55), breaks = round(seq(-0.5,0.5,0.25), digits = 2)),
      dv == "Relative\nimportance" ~ scale_x_continuous(limits = c(-0.21,0.65), breaks = seq(-0.2,0.6, 0.2) )
    )
  ) +
  
  
  geom_vline(xintercept = 0,
             lty = 2,
             colour = gray(1 / 2)) +
  ylab(element_text("party")) +
  theme(legend.title=element_blank()) +
  scale_fill_discrete(breaks=c("Conventional", "Bias-corrected", "Robust")) +
  scale_colour_manual(name="",
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  theme(legend.position="bottom") +
  ylab("Party") +
  # Reverse order of legend
  guides(color = guide_legend(reverse=TRUE),
         shape = guide_legend(reverse=TRUE)
  ) + 
  theme(panel.grid.minor.x = element_blank()) +
  
  theme(strip.text.x = element_text(size = 11, face = "bold"))


ggsave("plots/figure_g2.1a.pdf", width = 2100, height = 900, units = "px")



### Panel B -----------------------------------------------------------------

# Only results for Greens and AfD
results_party %>% 
  filter(party %in% c("green", "afd") ) %>% 
  mutate(party = factor(party, levels = c("green", "afd"), labels = c("Greens", "AfD"))) %>% 
  
  ggplot(aes(LATE, party)) +
  theme_bw() +
  geom_pointrange(
    aes(xmin = lower,
        xmax = upper,
        color = specification,
        shape = specification
    ),
    position = position_dodge(0.5),
    lwd = 0.6
  ) +
  scale_y_discrete(limits=rev) +
  facet_grid(cols = vars(dv), scales = "free_x") +
  ggh4x::facetted_pos_scales(
    x = list(
      dv == "Should represent own\nvoters in constituency" ~ scale_x_continuous(limits = c(-3.1,2.2), breaks = round(seq(-3,2,1), digits = 2)),
      dv == "Should represent all\ncitizens in constituency" ~ scale_x_continuous(limits = c(-3.1,2.2), breaks = round(seq(-3,2,1), digits = 2)),
      dv == "Relative\nimportance" ~ scale_x_continuous(limits = c(-0.51,2.55), breaks = seq(-0.5,2.5,0.5) )
    )
  ) +
  
  
  geom_vline(xintercept = 0,
             lty = 2,
             colour = gray(1 / 2)) +
  ylab(element_text("party")) +
  theme(legend.title=element_blank()) +
  scale_fill_discrete(breaks=c("Conventional", "Bias-corrected", "Robust")) +
  scale_colour_manual(name="",
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  theme(legend.position="bottom") +
  ylab("Party") +
  # Reverse order of legend
  guides(color = guide_legend(reverse=TRUE),
         shape = guide_legend(reverse=TRUE)
  ) + 
  theme(panel.grid.minor.x = element_blank()) +
  
  theme(strip.text.x = element_text(size = 11, face = "bold"))

ggsave("plots/figure_g2.1b.pdf", width = 2100, height = 900, units = "px")





## Figure G3.1: RDD estimates by election year -----------------------------

# NOTE: The code here follows the same structure as for Figure G2.1, with the only difference being that here it is subsetted to different election years instead of party affiliations

years <- c(2013,2017,2021)

specification <- c("Conventional", "Bias-corrected", "Robust")

results_own_years <- results_all_years <- results_relative_years <- array(NA, dim = c(length(years), 5, length(specification)))

colnames(results_own_years) <-
  colnames(results_all_years) <-
  colnames(results_relative_years) <-
  c("year",
    "LATE",
    "lower",
    "upper",
    "specification"
  )




for (i in 1:length(years)) {
  
  
  data_year <- post_data %>% filter(year == years[i])
  
  
  rdd_own_year <- rdrobust::rdrobust(y = data_year$rep_voters,
                                     x = data_year$margin,
                                     c = 0,
                                     p = 1,
                                     q = 2,
                                     covs = cbind(data_year$party_name),
                                     vce = "nn",
                                     cluster = data_year$cluster,
                                     all = TRUE
  )
  
  rdd_all_year <- rdrobust::rdrobust(y = data_year$rep_all,
                                     x = data_year$margin,
                                     c = 0,
                                     p = 1,
                                     q = 2,
                                     covs = cbind(data_year$party_name),
                                     vce = "nn",
                                     cluster = data_year$cluster,
                                     all = TRUE
  )
  
  rdd_relative_year <- rdrobust::rdrobust(y = data_year$relative_importance,
                                          x = data_year$margin,
                                          c = 0,
                                          p = 1,
                                          q = 2,
                                          covs = cbind(data_year$party_name),
                                          vce = "nn",
                                          cluster = data_year$cluster,
                                          all = TRUE
  )
  
  
  
  
  
  for (j in seq_along(specification)) {
    
    
    temp <- dplyr::if_else(j == 1, "tau.us", "tau.bc")
    
    results_own_years[i,1,j] <- years[i]
    results_own_years[i,2,j] <- rdd_own_year$Estimate[,temp]
    results_own_years[i,3,j] <- rdd_own_year$ci[j,1]
    results_own_years[i,4,j] <- rdd_own_year$ci[j,2]
    results_own_years[i,5,j] <- specification[j]
    
    
    results_all_years[i,1,j] <- years[i]
    results_all_years[i,2,j] <- rdd_all_year$Estimate[,temp]
    results_all_years[i,3,j] <- rdd_all_year$ci[j,1]
    results_all_years[i,4,j] <- rdd_all_year$ci[j,2]
    results_all_years[i,5,j] <- specification[j]
    
    results_relative_years[i,1,j] <- years[i]
    results_relative_years[i,2,j] <- rdd_relative_year$Estimate[,temp]
    results_relative_years[i,3,j] <- rdd_relative_year$ci[j,1]
    results_relative_years[i,4,j] <- rdd_relative_year$ci[j,2]
    results_relative_years[i,5,j] <- specification[j]
    
  }
  
}






temp_own <- rbind(results_own_years[,,1],
                  results_own_years[,,2],
                  results_own_years[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-("specification")), as.numeric) %>% 
  mutate(dv = "Own voters in constituency")

temp_all <- rbind(results_all_years[,,1],
                  results_all_years[,,2],
                  results_all_years[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-("specification")), as.numeric) %>% 
  mutate(dv = "All citizens in constituency")

temp_relative <- rbind(results_relative_years[,,1],
                       results_relative_years[,,2],
                       results_relative_years[,,3]) %>% 
  as.data.frame() %>% 
  mutate_at(vars(-("specification")), as.numeric) %>% 
  mutate(dv = "Relative importance")


results_years <- rbind(temp_own, temp_all, temp_relative) %>% 
  mutate(specification = factor(specification,
                                c("Robust", "Bias-corrected", "Conventional"),
                                labels = c("Robust", "Bias-corrected", "Conventional")
  )) %>%
  mutate(dv = factor(dv,
                     c("Own voters in constituency", "All citizens in constituency", "Relative importance"),
                     labels = c("Should represent own\nvoters in constituency", "Should represent all\ncitizens in constituency", "Relative\nimportance")
  )) %>% 
  mutate(year = factor(year, levels = c("2013","2017","2021")
  ))



results_years %>% 
  ggplot(aes(LATE, year)) +
  theme_bw() +
  geom_pointrange(
    aes(xmin = lower,
        xmax = upper,
        color = specification,
        shape = specification
    ),
    position = position_dodge(0.5),
    lwd = 0.6
  ) +
  scale_y_discrete(limits=rev) +
  facet_grid(cols = vars(dv), scales = "free_x") +
  ggh4x::facetted_pos_scales(
    x = list(
      dv == "Should represent own\nvoters in constituency" ~ scale_x_continuous(limits = c(-0.55,0.55), breaks = round(seq(-0.5,0.5,0.25), digits = 2)),
      dv == "Should represent all\ncitizens in constituency" ~ scale_x_continuous(limits = c(-0.55,0.55), breaks = round(seq(-0.5,0.5,0.25), digits = 2)),
      dv == "Relative\nimportance" ~ scale_x_continuous(limits = c(-0.21,0.81), breaks = seq(-0.2,0.8, 0.2) )
    )
  ) +
  geom_vline(xintercept = 0,
             lty = 2,
             colour = gray(1 / 2)) +
  ylab(element_text("Year")) +
  theme(legend.title=element_blank()) +
  scale_fill_discrete(breaks=c("Conventional", "Bias-corrected", "Robust")) +
  scale_colour_manual(name="",
                      values = c("Conventional" = viridis(3)[1], "Bias-corrected" = viridis(3)[2], "Robust" = viridis(3)[3])) +
  scale_shape_manual(name="",
                     values = c("Conventional" = 16, "Bias-corrected" = 2, "Robust" = 15)) +
  theme(legend.position="bottom") +
  # Reverse order of legend
  guides(color = guide_legend(reverse=TRUE),
         shape = guide_legend(reverse=TRUE)
  ) + 
  theme(panel.grid.minor.x = element_blank()) +
  
  theme(strip.text.x = element_text(size = 11, face = "bold"))


ggsave("plots/figure_g3.1.pdf", width = 1950, height = 1100, units = "px")

